1. Field of the Invention
The present invention relates generally to a method for subpixel registration of images and, more particularly, to a method for subpixel registration of satellite images of the earth""s surface to precisely determine environmental changes of the surface.
2. Prior Art
Satellite images of the earth""s surface are used to study environmental changes by comparing images of particular areas taken over a period of time. In order to measure the changes due to environmental factors, it is necessary to register the images to each other so that a region in one image is aligned perfectly with the corresponding view of the same region in a second image. When aligned properly, the pixels in the first image lie directly on top of the corresponding pixels in the second image. Image registration may involve translation (side-to-side or up-down movement), rotation, re-scaling by uniform shrinking or enlarging, and more complicated warping of the images.
There are those in the art who have shown that registration has to be very accurate in order to measure environmental change because the disparity observed due to mis-registration cannot be distinguished from changes due to the environment. Therefore, sensitivity to mis-registration is very large, which places a requirement on registration that it be done to subpixel precision.
It is known in the art that the cross-power spectrum of two band-limited images is the Fourier transform of the function sin(Πx) sin(Πy)/Π2x y, where x and y are the translational offset of one image with respect to the other. When the cross power spectrum is inverted, one can use two or three samples near the peak and solve for the values of x and y. However, satellite images do not satisfy the band-limiting constraint required for this process. Band-limited images are necessarily very smooth, and do not contain sharp features. Satellite images tend to have many sharp features and contain a good deal of energy in the high frequencies. Thus this process yields inaccurate results when applied to such images.
One method of the prior art proposes to fit a parabola to the correlation peak, and use this fit to obtain an estimate of the peak position. As an alternative it also proposes to fit a Gaussian to the correlation peak. This prior art method also describes a method that takes a somewhat different tack. This scheme computes the phase differences between the Fourier transforms of the two images for the four integer-valued translation amounts nearest to the estimate of the true translation. The scheme computes the final estimate of the true translation by finding the linear interpolation of the four phase-relations that best approximates the observed phase relation.
Others have reported variations of this idea. It is well-known that the phase relation of the Fourier transforms of two images that are slightly displaced with respect to each other contain information about the displacement. Aliasing effects on the transforms make it difficult to extract that information. There are those in the art who have proposed to use the low frequency portion of the spectrum, because that portion is less likely to be aliased.
Another portion of prior art is in the area of time delay estimation in radar and sonar returns in which a maximum likelihood estimator for signal delay is developed, which corresponds to the problem to image translation. This method assumes that both noise and signal are band-limited, and this allows the signal to be interpolated accurately between sample points. In the presence of noise, the registrations are generally not perfect for any delay time, and the problem becomes one of choosing a delay time that is the most likely one in the sense that the noise that explains the discrepancies is the most likely to have been observed. This process uses recursive demodulation to register two signals. The recursion alternates between estimating the delay and estimating the signal. It uses a maximum likelihood criterion that depends on the noise present in the observations. A related scheme for dealing with correlated noise that involves the calculation of the covariance matrix to remove the effects of noise. A later described method produces maximum likelihood estimates from a scheme that is an alternative to the above discussed expansion, which in turn is based on the covariance matrix of the samples. For narrow bandwidth signals, it has been described to use an adaptive delay algorithm that uses a steepest descent technique to find a delay that produces a minimum mean square error.
Others in the art approached this problem as one in adaptive filtering. They describe modifying a filter in real time to adapt to an incoming signal whose delay changes dynamically. They also describe the static and dynamic behavior, respectively, of their tracking scheme.
Phase estimates in radar and sonar systems play a role analogous to estimating time delays of signals. The phase relation of the Fourier transforms of time-shifted signals or pixel-shifted images reveals the amount of the shift. Those in the art described a means for estimating phase by first derivative of the phase, and then performing an integral of the derivative function. This process eliminates ambiguity of multiples of 2Π that occur when phase is estimated directly. Others presented an alternative method for estimating phase that is done in the Fourier domain on a symmetrized estimate of the phase.
In the area of image registration, still others used a multi-resolution approach to subpixel registration of images based on image features. Their method allows image warping to be one of the image transformations by which two images differ. Because warping is permitted, an ideal registration always exists, and the problem becomes one of deciding which of the possible transforms that give ideal results is the most likely to explain the difference in the images. In some applications, particularly in digital video applications, more than two images of the same scene are available with the scenes differing by small movements. There are those in the art that describe methods for using multiple images to construct a single super-resolution image from which displacements can be calculated to subpixel accuracy. These schemes are based on using averages and interpolations in the pixel domain to construct the super-resolution image.
A method is also known in the art that enables one to predict in advance what the accuracy of subpixel registration can be for a given image. Their scheme is based on the analysis of features within an image. Their registration method matches features, particularly line segments, to find the best match, and they simultaneously give an estimate of the precision of the registration.
The problem of discovering a displacement vector field that describes how two images differ is a registration problem that arises when the differences can be attributed to complex warping transformations. Instead of explaining the differences by a single translation or rotation, or a combination of the two, the differences must be explained by motion directions that could be different at every pixel in a scene. There are those in the art who describe an algorithm based on expectation maximization and have developed a method for this problem when the illumination source is quantum limited. They too use a maximum likelihood estimator to solve this problem.
Therefore it is an object of the present invention to provide a method for subpixel registration of images which yields accurate results when applied to satellite images.
It is yet a further object of the present invention to provide a method for subpixel registration of images which are not based on using averages and interpolations in the pixel domain.
Accordingly, a first embodiment of the methods of the present invention for registering first and second images which are offset by an x and/or y displacement in sub-pixel locations is provided. The method of the first embodiment comprising the steps of: multiplying the first image by a window function to create a first windowed image; transforming the first windowed image with a Fourier transform to create a first image Fourier transform; multiplying the second image by the window function to create a second windowed image; transforming the second windowed image with a Fourier transform to create a second image Fourier transform; computing a collection of coordinate pairs from the first and second image Fourier transforms such that at each coordinate pair the values of the first and second image Fourier transforms are likely to have very little aliasing noise; computing an estimate of a linear Fourier phase relation between the first and second image Fourier transforms using the Fourier phases of the first and second image Fourier transforms at the coordinate pairs in a minimum-least squares sense; and computing the displacements in the x and/or y directions from the linear Fourier phase relationship.
Also provided is a variation of the method of the first embodiment of the present invention for registering first and second images which are offset by an x and/or y displacement in sub-pixel locations. The method comprising the steps of: registering the first and second images to the nearest pixel location; and registering the first and second images to the nearest sub-pixel location, wherein the registering of the first and second images to the nearest sub-pixel location comprises the steps of: multiplying the first image by a window function to create a first windowed image; transforming the first windowed image with a Fourier transform to create a first image Fourier transform; multiplying the second image by the window function to create a second windowed image; transforming the second windowed image with a Fourier transform to create a second image Fourier transform; computing a collection of coordinate pairs from the first and second image Fourier transforms such that at each coordinate pair the values of the first and second image Fourier transforms are likely to have very little aliasing noise; computing an estimate of a linear Fourier phase relation between the first and second image Fourier transforms using the Fourier phases of the first and second image Fourier transforms at the coordinate pairs in a minimum-least squares sense; and computing the displacements in the x and/or y directions from the linear Fourier phase relationship.
In a preferred implementation of the methods of the first embodiment of the present invention, the step of computing a collection of coordinate pairs comprises a conjunction of two or more of the following: fixing the coordinate pairs to lie within a fixed distance of a central peak; fixing the Fourier magnitude of the first image Fourier transform and of the second image Fourier transform to be among the highest predetermined percent of the magnitudes within a fixed distance of a central peak; and setting the ratio of the Fourier magnitude of the first image transform to the Fourier magnitude of the second image transform to lie between 1/(1+xcfx81) and 1+xcfx81, where xcfx81 is a positive number.
Also provided is a computer program product embodied in a computer-readable medium for registering first and second images which are offset by an x and/or y displacement in sub-pixel locations. The computer program product comprises: computer readable program code means for multiplying the first image by a window function to create a first windowed image; computer readable program code means for transforming the first windowed image with a Fourier transform to create a first image Fourier transform; computer readable program code means for multiplying the second image by the window function to create a second windowed image; computer readable program code means for transforming the second windowed image with a Fourier transform to create a second image Fourier transform; computer readable program code means for computing a collection of coordinate pairs from the first and second image Fourier transforms such that at each coordinate pair the values of the first and second image Fourier transforms are likely to have very little aliasing noise; computer readable program code means for computing an estimate of a linear Fourier phase relation between the first and second image Fourier transforms using the Fourier phases of the first and second image Fourier transforms at the coordinate pairs in a minimum-least squares sense; and computer readable program code means for computing the displacements in the x and/or y directions from the linear Fourier phase relationship.
Yet still provided is a program storage device readable by machine, tangibly embodying a program of instructions executable by the machine to perform method steps for registering first and second images which are offset by an x and/or y displacement in sub-pixel locations. Wherein the method comprises the steps of: multiplying the first image by a window function to create a first windowed image; transforming the first windowed image with a Fourier transform to create a first image Fourier transform; multiplying the second image by the window function to create a second windowed image; transforming the second windowed image with a Fourier transform to create a second image Fourier transform; computing a collection of coordinate pairs from the first and second image Fourier transforms such that at each coordinate pair the values of the first and second image Fourier transforms are likely to have very little aliasing noise; computing an estimate of a linear Fourier phase relation between the first and second image Fourier transforms using the Fourier phases of the first and second image Fourier transforms at the coordinate pairs in a minimum-least squares sense; and computing the displacements in the x and/or y directions from the linear Fourier phase relationship.
Still yet provided is a second embodiment of the methods of the present invention for registering first and second images which are offset by an x and/or y displacement in sub-pixel locations. The second embodiment of the method comprises the steps of: (a) Fourier transforming the first image to produce a first image Fourier transform; (b) Fourier transforming the second image to produce a second image Fourier transform; (c) summing the energy of the first image Fourier transform and the energy of the second image Fourier transform and dividing the sum by 2 to produce an average image energy; (d) selecting one or more first displacements for an iterative search; (e) using each first displacement, the first-image Fourier transform, and the second-image Fourier transform to compute a first baseband-plus-aliasing Fourier transform; (f) comparing the energy of the first baseband-plus-aliasing Fourier transform to the average image energy; (g) stopping the search if the energy comparisons are within a fixed precision; (h) continuing the search if the energy comparisons are not within the fixed precision by selecting a new displacement for the search; (i) computing a second baseband-plus-aliasing Fourier transform using the new displacement, first-image Fourier transform, and the second-image Fourier transform; (j) repeating steps (h) and (i) until the fixed precision is reached or a criterion for stopping the search unsuccessfully is satisfied.